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ABSTRACT 

A residual planetesimal disk of mass 10-100M® remained in the outer solar system following the 
birth of the giant planets, as implied by the existence of the Oort cloud, coagulation requirements 
for Pluto, and inefficiencies in planet formation. Upon gravitationally scattering planetesimal debris, 
planets migrate. Orbital migration can lead to resonance capture, as evidenced here in the Kuiper 
and asteroid belts, and abroad in extra-solar systems. Finite sizes of planetesimals render migration 
stochastic ("noisy"). At fixed disk mass, larger (fewer) planetesimals generate more noise. Extreme 
noise defeats resonance capture. We employ order-of-magnitude physics to construct an analytic 
theory for how a planet's orbital semi-major axis fluctuates in response to random planetesimal scat- 
terings. The degree of stochasticity depends not only on the sizes of planetesimals, but also on their 
orbital elements. We identify the conditions under which the planet's migration is maximally noisy. 
To retain a body in resonance, the planet's semi-major axis must not random walk a distance greater 
than the resonant libration width. We translate this criterion into an analytic formula for the reten- 
tion efficiency of the resonance as a function of system parameters, including planetesimal size. We 
verify our results with tailored numerical simulations. Application of our theory reveals that capture 
of Resonant Kuiper belt objects by a migrating Neptune remains effective if the bulk of the primordial 
disk was locked in bodies having sizes < 0(100) km and if the fraction of disk mass in objects with 
sizes > 1000 km was less than a few percent. Coagulation simulations produce a size distribution of 
primordial planetesimals that easily satisfies these constraints. We conclude that stochasticity did not 
interfere with nor modify in any substantive way Neptune's ability to capture and retain Resonant 
Kuiper belt objects during its migration. 

Subject headings: celestial mechanics — Kuiper belt — diffusion — planets and satellites: formation — 
solar system: formation 



1. INTRODUCTION 

Planet formation by coagulation of planetesimals is 
not perfectly efficient — it leaves behind a residual disk of 
solids. Upon their coalescence, the outer planets of our 
solar system were likely embedded in a 10-100M® disk of 
rock and ice containing the precursors of the Oort cloud 
(Doncs et al. 2004) and the Kuiper belt (see the reviews 
by Chiang et al. 2006; Cruikshank et al. 2006; Levison et 
al. 2006). The gravitational back-reaction felt by planets 
as they scatter and scour planetesimals causes the plan- 
ets to migrate (Fernandez & Ip 1984; Murray et al. 1998; 
Hahn & Malhotra 1999; Gomes, Morbidelli, & Levison 
2004). Neptune is thought to have migrated outward 
and thereby trapped Kuiper belt objects (KBOs) into 
its exterior mean-motion resonances, both of low-order 
such as the 3:2 (Malhotra 1995) and of high-order such 
as the 5:2 (Chiang et al. 2003; Hahn & Malhotra 2005). 
Likewise, Jupiter's inward migration may explain the ex- 
istence of Hilda asteroids in 2:3 resonance with the gas gi- 
ant (Franklin et al. 2004). A few pairs of extra-solar plan- 
ets, locked today in 2:1 resonance (Vogt et al. 2005; Lee et 
al. 2006), may have migrated to their current locations 
within parent disks composed of gas and/or planetesi- 
mals. Orbital migration and resonant trapping of dust 
grains may also be required to explain non-axisymmetric 
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structures observed in debris disks surrounding stars 10- 
100 Myr old (e.g., Wyatt 2003; Meyer et al. 2006). 

Only when orbital migration is sufficiently smooth and 
slow can resonances trap bodies. The slowness crite- 
rion requires migration to be adiabatic: Over the time 
the planet takes to migrate across the width of the res- 
onance, its resonant partner must complete at least a 
few librations. Otherwise the bodies speed past reso- 
nance (e.g., Dermott, Malhotra, & Murray 1988; Chiang 
2003; Quillen 2006). Smoothness requires that changes 
in the planet's orbit which are incoherent over timescales 
shorter than the libration time do not accumulate un- 
duly. Orbital migration driven by gravitational scatter- 
ing of discrete planetesimals is intrinsically not smooth. 
A longstanding concern has been whether Neptune's mi- 
gration was too "noisy" to permit resonance capture and 
retention (see, e.g., Morbidelli, Brown, & Levison 2003). 
In N-body simulations of migration within planetesimal 
disks (Hahn & Malhotra 1999; Gomes et al. 2004; Tsiga- 
nis et al. 2005), N ~ O(10 4 ) is still too small to produce 
the large, order-unity capture efficiencies seemingly de- 
manded by the current census of Resonant KBOs. 

At the same time, the impediment against resonance 
capture introduced by inherent stochasticity has been ex- 
ploited to explain certain puzzling features of the Kuiper 
belt, most notably the Classical (non-Resonant) belt's 
outer truncation radius, assumed to lie at a heliocentric 
distance of ~48 AU (Trujillo & Brown 2001; Levison & 
Morbidelli 2003). If Neptune's 2:1 resonance captured 
KBOs and released them en route, Classical KBOs could 
have been transported ("combed") outwards to popu- 
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late the space interior to the final position of the 2:1 
resonance, at a semi-major axis of 47.8 AU (Levison & 
Morbidelli 2003). As originally envisioned, this scenario 
requires that ~3M^ be trapped inside the 2:1 resonance 
so that an attendant secular resonance suppresses growth 
of eccentricity during transport. It further requires that 
the degree of stochasticity be such that the migration 
is neither too smooth nor too noisy. Whether these re- 
quirements were actually met remain open questions. 3 

Stochastic migration has also been studied in gas disks, 
in which noise is driven by density fluctuations in tur- 
bulent gas. Laughlin, Steinacker, & Adams (2004) and 
Nelson (2005) propose that stochasticity arising from 
gas that is unstable to the magneto-rotational instabil- 
ity (MRI) can significantly prolong a planet's survival 
time against accretion onto the parent star. The spec- 
trum of density fluctuations is computed by numerical 
simulations of assumed turbulent gas. 

In this work, we study stochastic changes to a planet's 
orbit due to planetesimal scatterings. The planet's 
Brownian motion arises from both Poisson variations 
in the rate at which a planet encounters planetcsimals, 
and from random fluctuations in the mix of planet- 
planetesimal encounter geometries. How does the vigor 
of a planet's random walk depend on the masses and or- 
bital properties of surrounding planetesimals? We an- 
swer this question in j}2l by constructing an analytic 
theory for how a migrating planet's semi-major axis 
fluctuates about its mean value. We employ order-of- 
magnitude physics, verifying our assertions whenever fea- 
sible by tailored numerical integrations. Because the 
properties of planetesimal disks during the era of plane- 
tary migration are so uncertain, we consider a wide vari- 
ety of possibilities for how planetesimal semi-major axes 
and eccentricities are distributed. One of the fruits of 
our labors will be identification of the conditions under 
which a planet's migration is maximally stochastic. 

Apportioning a fixed disk mass to fewer, larger plan- 
etesimals renders migration more noisy. How noisy is too 
noisy for resonance capture? What limits can we place 
on the sizes of planetesimals that would keep capture 
of Resonant KBOs by a migrating Neptune a viable hy- 
pothesis? These questions are answered in where we 
write down a simple analytic formula for the retention 
efficiency of a resonance as a function of disk proper- 
ties, including planetesimal size. Quantifying the size 
spectrum of planetesimals is crucial for deciphering the 
history of planetary systems. Many scenarios for the evo- 
lution of the Kuiper belt implicitly assume that most of 
the mass of the primordial outer solar system was locked 
in planetesimals having sizes of 0(100) km, like those ob- 
served today (see, e.g., Chiang et al. 2006 for a critique 
of these scenarios). By contrast, coagulation simulations 
place the bulk of the mass in bodies having sizes of 0(1) 
km (Kenyon & Luu 1999). For ice giant formation to 
proceed in situ in a timely manner in the outer solar sys- 

3 While Classical KBOs do have semi-major axes that extend 
up to 48 AU, the distribution of their perihelion distances cuts off 
sharply at distances closer to 45 AU (see, e.g, Figure 2 of Chiang et 
al. 2006). Interpreted naively (i.e., without statistics), the absence 
of bodies having perihelion distances of 45-48 AU and eccentricities 
less than ~0.1 smacks of observational bias and motivates us to re- 
visit the problem of whether an edge actually exists, or at least 
whether the edge bears any relation to the 2:1 resonance. 



tern, most of the primordial disk may have to reside in 
small, sub- km bodies (Goldreich, Lithwick, & Sari 2004). 

In 21 m addition to summarizing our findings, we ex- 
tend them in a few directions. The main thrust of this 
paper is to analyze how numerous, small perturbations 
to a planet's orbit accumulate. We extend our analysis 
in 21 1° quantify the circumstances under which a sin- 
gle kick to the planet from an extremely large planetesi- 
mal can disrupt the resonance. We also examine pertur- 
bations exerted directly on Resonant KBOs by ambient 
planetesimals. 

2. STOCHASTIC MIGRATION: AN 
ORDER-OF-MAGNITUDE THEORY 

We assume the planet's eccentricity is negligibly small. 
We decompose the rate of change of the planet's semi- 
major axis, d p , into average and random components, 

d p — ^p,avg "F d P; rnd • (1) 

The average component ( "signal" ) arises from any global 
asymmetry in the way a planet scatters planetesimals, 
e.g., an asymmetry due to systematic differences between 
planetesimals inside and outside a planet's orbit. The 
random component ( "noise" ) results from chance varia- 
tions in the numbers and orbital elements of planetesi- 
mals interacting with the planet. By definition, d Pjrn d 
time-averages to zero. We assume that a P:av Jt) is a 
known function of time t, and devote all of J21 to the 
derivation of d p , rn d- 

Each close encounter between the planet and a single 
planetesimal lasting time At c causes the planet's semi- 
major axis to change by Aet p . Expressions for Aa p and 
Af e depend on the planetesimal's orbital elements. We 
define x = a — a p as the difference between the semi- 
major axes of the planetesimal and of the planet, b > 
as the impact parameter of the encounter, and u ~ eQa 
as the planetesimal's random (epicyclic) velocity, where 
a, e, and f2 are the semi-major axis, eccentricity, and 
mean angular velocity of the planetesimal, respectively. 
We assume that |x| < a p . Encounters unfold differently 
according to how \x\ and b compare with the planet's Hill 
radius, 

and according to how u compares with the Hill velocity, 

v H = f> p i? H • (3) 

Here M p and M* are the masses of the planet and of 
the star, respectively, and S7 p is the angular velocity of 
the planet. See Table IaT1 for a listing of frequently used 
symbols. 

In t l2.ll we calculate Aa p for a single encounter with 
a planetesimal having |x| > i?H- In H2.2I we repeat the 
calculation for |x| < i?H- In M2.3I we provide formulae 
for the root-mean-squared (RMS) random velocity due 
to cumulative encounters, (dp m( i) ' and identify which 
cases of those treated in S ^2.1hl2~2l potentially yield the 
strongest degree of stochasticity in the planet's migra- 
tion. 
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2.1. Single Encounters with |x| > i?u- Non-Horseshoes 

We calculate the change in the planet's semi-major 
axis, Aa p , resulting from an encounter with a single plan- 
etesimal having |x| > i?H- We treat planetesimals on 
orbits that do not cross that of the planet in M2.1.1I and 
those that do cross in M2.1.2I Throughout, A refers to the 
change in a quantity over a single encounter, evaluated 
between times well before and well after the encounter. 

2.1.1. N on- Crossing Orbits 

Planetesimals on orbits that do not cross that of the 
planet have 

|x| > ae , (4) 

which corresponds to 

\x\/Ru > u/v H . (5) 

Our plan is to relate Aa p to Ax by conservation of en- 
ergy, calculate Ae using the impulse approximation, and 
finally generate Ax from Ae by conservation of the Ja- 
cobi integral. 
By conservation of energy, 



A 



GM*M D GM*m 



2a r 



2a 







(6) 



where m <C M p is the mass of the planetesimal. We have 
dropped terms that account for the potential energies of 
the planet and of the planetesimal in the gravitational 
field of the ambient disk. These are small because the 
disk mass is of order M p -C M* and because the disk 
does not act as a point mass but is spatially distributed. 
Equation JSJ implies 



Aa c 



m /Op\ 
M D \a J 



Aa . 



Since |Aa p | <C |Aa| and a p 



Aa r 



~ a, we have Ax < 
m 

M 



■Ax . 



(7) 

Aa and 

(8) 



The impulse imparted by the planet changes the ec- 
centricity of the planetesimal by Ae. An encounter for 
which | a; | is more than a few times i?H imparts an impulse 



per mass 



An ~ ±^At c . (9) 



The impact parameter b is limited by 

\x\-ae<b<\x\+ae . (10) 

Since ae < |x|, 

b~\x\. (11) 

Because the relative speed due to Keplerian shear, 
(3/2)Q p |x|, is larger than u, the relative speed during 
encounter is dominated by the former, and 

2b _ 4 1 
(3/2)ft p 6 ~~ 30p" ~ Op" ' 



AU 



(12) 



Since At e is about one-fifth of an orbital period, the im- 
pulse approximation embodied in @ should yield good 

4 The impulse to the planetesimal changes both u and the plan- 
etesimal's Keplerian shearing velocity, — (3/2)fia;. In the non- 
crossing case, |Au| > [A(f2a;)[. 



order-of-magnitude results. The change in the eccentric- 
ity of the planetesimal is hence 



Ae 



At 



M* V x J 



(13) 



When |Ae| < (the pre-encounter) e, the change Ae can 
be either positive or negative, depending on the true 
anomaly of the planetesimal at the time of encounter. 
If |Ae| > e, then Ae > 0. When \x\ ~ i? H , |Ae| attains 
its maximum value of ^Mp/M*) 1 / 3 ; i.e., Am ~ «h- 5 

To calculate the corresponding change in the planetes- 
imal's semi-major axis, Ax, we exploit conservation of 
the Jacobi integral, Cj. That a conserved integral ex- 
ists relies on the assumption that in the frame rotating 
with the planet, the potential (having a centrifugal term 
plus gravitational contributions due to the star, planet, 
and disk) is time-stationary; the Jacobi integral is simply 
the energy of the planetesimal (test particle) evaluated 
in that frame. To the same approximation embodied in 
Equation JHJl, 



1 



Cj=E- ftp J 
GNU 



2{a/a p ) 



+ V(a/o P )(l-e 2 ) 



(14) 



far from encounter, where E and J are the energy and 
angular momentum per mass of the planetesimal, respec- 
tively. Taking the differential of (|14fl yields, to leading 
order, 



3 /Ax 

4 V a„ 



3 x 
2 a n 



1 9 \ Ax 
-e 



A(e 2 ) = 0. (15) 



Since |x|/a p > e > e 2 (non-crossing condition) and 
A(e 2 ) < (x/ap) 2 (by Equation an d the condition 
|x| > -Rh), Equation l(TK|l reduces to 



2a 2 

Ax~-^A(e 2 ) . 
3x 

We combine Equations © and i|16[l to find 
Aa p ~ - — ^A e 2 . 

M D X 



(16) 



(17) 



Equation l|17|) takes two forms depending on how |Ae| 
compares with (the pre-encounter) e. If 



-""It 1 ' 2 



(18) 



then |Ae| < e, A(e 2 ) ~ 2eAe, and Equation JTJl be- 
comes 



a ma 
Aa D ~ =F re . 

p M* x 3 



(19) 



5 When |x| < 2_Rjj, the encounter pulls the planetesimal into the 
planet's Hill sphere. The planetesimal accelerates in a complicated 
way and exits the Hill sphere in a random direction with u of 
order the planet's escape velocity at the Hill radius, ujj (Petit & 
Henon 1986). The planetesimal's eccentricity is boosted by Ae ~ 
(Mp/M,) 1 / 3 . The encounter time is typically the time required to 
complete a few orbits around the planet, At c ~ 2n/ Q . Since Af c 
and Ae match, to order of magnitude, Equations 1121 and 1131 
for \x\ ~ -Rh, we do not treat _Rh < \x\ < 2/?h as an explicitly 
different case. 
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The right-hand side is extremized for \x\ ~ R^vn/u) 1 ^ 2 : 

/ \ 5 / 2 
to „ / a n e \ to „ , 

max |Aa p | ~ — i? H [ ] < — -Rh , (20) 



Mn 



valid for e < Rn/a p (non-crossing). 
On the other hand, if 



i? H < x < Rn (— ) 



(21) 



then |Ae| > e and A(e 2 ) ~ (Ae) 2 , and Equation itTTll 
becomes 



Aa n 



mM p a p 
' M? x$ 



(22) 



Equation 122|) agrees with the more careful solution of 
Hill's problem by Henon & Petit (1986). The right-hand 
side is extremized for Ixl ~ Rr: 



max \Aa 



TO 



(23) 



In summary, if x/ Rr > u/vn, then (a) the planetesi- 
mal's orbit does not cross ("nc") that of the planet, (b) 



At c 



fit 



(24) 



and (c) 



Aa 



p,ncl 



a ma 



" M 2 a?' 

if i?H<x<(« H /u) 1/2 i? H ; 

4 



if * > (W«) 1/2 i?H- 



(25) 



2.1.2. Crossing Orbits 



Encounters with planetesimals on orbits that cross that 
of the planet, i.e., those with 



\x\/Rh < u/v K , 



(26) 



differ from encounters with non-crossing planetesimals 
in two key respects. First, the relative velocity of the 
two bodies is dominated by the planetesimal's random 
(epicyclic) velocity rather than the Keplerian shear. Sec- 
ond, the planetesimal's impact parameter, b, may differ 
significantly from |x|. The impact parameter may take 
any value 



< b < ae , 



(27) 



where o m j n is the impact parameter below which the plan- 
etesimal collides with the planet. Because crossing orbits 
allow for encounters with many different geometries, out- 
comes of these encounters can vary dramatically. Here we 
restrict ourselves to estimating the maximum |Aa p J that 
can result from an orbit-crossing encounter. In H2.3.2I 
we argue this restriction is sufficient for our purposes. 

When u > i>h, the eccentricity of the planetesimal 
can change by at most |Ae| ~ e. Such a change corre- 
sponds to an order-unity rotation of the direction of the 
planetesimal's random velocity vector, and requires that 



b < GM p /u 2 . The change in the planetesimal's specific 
energy over the encounter is approximately 



A - 



GM* 
2a 



GAL 



(28) 



where v is the velocity of the planetesimal relative to 
the star (in an inertial frame of reference) and r is the 
distance between the planetesimal and the star. Now Ar 
for an encounter with b < GM p /u 2 is at most GM p /u 2 < 
i?H and A(v 2 ) is of order Vlau. Since 

GM P , . 

, (29) 



Oftii > n 2 aRa > n 2 c 



the second term on the right-hand side of Equation l|28|l 
is negligible compared to the first, and the maximum 
|Ao| over an encounter is 



max lAol 



GM, 



■Qau ~ ae 



By Equation J7J) and a ~ a 
max I Aa 



(30) 



(31) 



We have verified Equation l|31l) by numerical orbit in- 
tegrations. We could also have arrived at Equation l|31|l 
through Equation i|15|) . which yields |Ax| ~ a p e for cross- 
ing orbits when |Ae| ~ e. 

2.2. Single Encounters with \x\ < i?n-' Horseshoes 

When |x| < Rh, planetesimals can occupy horseshoe 
orbits. A planetesimal on a horseshoe orbit for which 
|at| ~ i?H encounters the planet on a timescale somewhat 
shorter than the orbital period; by the impulse approxi- 
mation, such a planetesimal kicks the planet such that 

(32) 



A m % 



where we have momentarily restricted consideration to 
planetesimals having sub- Hill eccentricities (e < Rn/a p ). 
From Henon and Petit (1986), 



3 x 2 ■ 



valid for x not too far below i?H- Then 
'A./„ 



m x 



(33) 



(34) 



M p Rl " 

The kick is maximal for maximum \x\ — R^: 

Tfi 

max|Aa p | — ' ^ 35 ^ 

This is the same maximum as was derived for the |x| ~ 
i?H, non-crossing case; see Equation (|23|l . Thus, a co- 
orbital ring of planetesimals on horseshoe orbits with 
sub-Hill eccentricities increases the stochasticity gener- 
ated by planetesimals on non-horseshoe, non-crossing or- 
bits by a factor of at most order unity (under the assump- 
tion that disk properties are roughly constant within sev- 
eral Hill radii of the planet). For this reason, and also be- 
cause the horseshoe region may well have been depleted 
of planetesimals compared to the rest of the disk, we 
omit consideration of co-orbital, sub-Hill planetesimals 
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for the remainder of the paper, confident that the error 
so incurred will be at most order unity. 

What about planctesimals on horseshoe orbits with 
super-Hill eccentricities (e > Rn/a p )7 Upon encounter- 
ing the planet, such objects can have their semi-major 
axes changed by |Aa| > i?n — whereupon they are ex- 
pelled from the 1:1 horseshoe resonance. Because highly 
eccentric, horseshoe resonators are unstable, we neglect 
consideration of them for the rest of our study. 

2.3. Multiple Encounters: Cumulative Stochasticity 

We now extend our analysis from individual encoun- 
ters to the cumulative stochasticity generated by a disk 
with surface density E m in planetesimals of a single mass 
m. Note that E m need not equal the total surface den- 
sity £ (integrated over all possible masses m). We will 
consider size distributions in ^3.61 We consider plan- 
etesimals with sub-Hill (u < «h) velocities ( W2. 3. 1(1 sep- 
arately from those with super-Hill (u > uh) velocities 
( H2.3.2JI . Sub-Hill (non- horseshoe) planetesimals always 
occupy non-crossing orbits. Super-Hill planetesimals can 
be crossing or non-crossing. 

2.3.1. Sub-Hill Velocities (u < vn) 

Consider planetesimals with sub-Hill velocities located 
a radial distance x away from the planet (|a:| > i?n)- 
Since u < «h, the speeds of planetesimals relative to 
the planet are determined principally by Keplerian shear 
(Equation ^2])> an d the scale height of the planetesi- 
mals is less than i?n- The planet encounters (undergoes 
conjunctions with) such planetesimals at a mean rate 



m 



n x 2 



(36) 



as is appropriate for encounters in a two-dimensional ge- 
ometry. Over a time interval At, the planet encounters 

N = NAt such planetesimals on average. Systematic 
trends in N with x — say, systematically more objects 
encountered interior to the planet's orbit than exterior 
to it — cause the planet to migrate along an average tra- 
jectory with velocity a p , a vg- 

Random fluctuations in (a) the number of planetesi- 
mals encountered per fixed time interval and (b) the mix 
of planetesimals' pre-encounter orbital elements cause 
the planet to random walk about this average trajectory. 
Contribution (a) is straightforward to model. The prob- 
ability that the planet encounters N objects located a 
distance x away in time At is given by Poisson statistics: 



P(N) = *le-» 



(37) 



The variance in N is 



a 2 N = {(N-N) 2 ) = N . (38) 

Fluctuations in N drive the planet either towards or away 
from the star with equal probability and with typical 
speed 



,1/2 



|Aa p |-i/2^ 



p,rnd/ '~ fa l (39) 

hereafter the root-mean-squared (RMS) speed. While 

cx 1/ V At, the distance random walked 
1 / 2 At cx \/Al. 



(a 2 W 2 

\ u p,rnd/ 



Our assumption of Poisson statistics is reasonable. In 
the sub-Hill planet-planetesimal encounter re- 

quires a time At e ~ l/f2 p (Equation |12|) to complete. 
Encounters separated by more than At e are uncorrelated 
with one another, at least until the planet completes one 
revolution with respect to the surrounding disk, i.e., at 
least until a synodic time t syn ~ 47ra p /(3f2 p |2;|) elapses. 
After a synodic period, it is possible, in principle, for the 
planet to essentially repeat the same sequence of encoun- 
ters that it underwent during the last synodic period. We 
assume in this paper that this does not happen — that 
the orbits of planetesimals interacting with the planet 
are randomized on a timescale t v< \ z < t syn . We ex- 
pect this inequality to be enforced by a combination of 
(i) randomization of planetesimal orbits due to encoun- 
ters with the planet (e.g., encounters within the chaotic 
zone of the planet [Wisdom 1980]), (ii) phase mixing of 
planetesimals due to Keplerian shear (which occurs on 
timescale t syn for planetesimals distributed between x 
and ^2x), (iii) gravitational interactions between plan- 
etesimals, and (iv) physical collisions between planetesi- 
mals. As long as t r( jz < ^syn, we are free to choose At to 
be anything longer than At c . 6 

Contribution (b) is difficult to model precisely since 
we do not know how orbital elements of planetesimals 
are distributed. These distributions are unlikely to be 
governed by simple Poisson or Gaussian statistics (see, 
e.g., Ida & Makino 1992; Rafikov 2003; Collins & Sari 
2006). Nonetheless, neglecting contribution (b) will not 
lead to serious error. Suppose the planetesimals' pre- 
encounter elements are distributed such that the frac- 
tional variation in each element is at most of order unity 
(e.g., the planetesimal eccentricities span a range from 
e/2 to 2e at most). Then the central limit theorem en- 
sures that the noise introduced by random sampling of 
orbital elements is at most comparable to the noise in- 
troduced by random fluctuations in the encounter rate. 
Consider, for example, noise that arises from random 
sampling of e in the case where Aa p = Aa PjnC 2 (Equation 

|25|). For an encounter rate fixed at N, the planet's semi- 
major axis a p changes over time interval At by N x Aa p , 

where Aa p is the mean of N = NAt sampled values of 
Aa p . If the dispersion in e for individual planetesimals 
is <j e and the mean eccentricity sampled over N values 
is e, then the dispersion in the sampled mean eccentric- 

1/2 

ity is <7e ~ a e /N by the central limit theorem. For 
Aa p = Aa P)nC 2 cx e, the dispersion in Aa p is \Aa p \o~e/e. 
The planet's RMS speed generated purely from random 
sampling of e is 



\l/2 



N IA I °* 

Ai |A ^' ~ 



\Aa 



At 



PI fTe 



(40) 



*p,rnd/ 



which is at most comparable to the RMS speed gener- 
ated purely from random sampling of N (Equation |39jl. 

6 If t r( j z > t syn , then At > t rdz and the right-hand side of 
Equation 1391 is multiplied by y 7 t r ^/t syn . The planet's motion 
is more stochastic in this case because over t r( j z , correlated inter- 
actions with planetesimals do not cancel each other as much as 
uncorrelated interactions would. Later, since we will be interested 
in stochastic perturbations to mean-motion resonant particles, we 
will require t T ^ z < tin,, where t^t, is the libration period within 
resonance. 
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as desired. Our original supposition that max(cr e /e) ~ 1 
leads to no important loss of generality; if the distri- 
bution in e were bi-modal, for example, we could treat 
each population separately and add the resultant RMS 
speeds in quadrature. Similar results obtain for random 
sampling of other elements such as x. For simplicity, 
we hereafter treat explicitly only fluctuations in the en- 
counter rate [(Equation knowing that the noise so 
calculated will be underestimated by a factor of at most 
order unity. 

Expression 139|) measures the contribution to the RMS 
speed from planetesimals located a distance x away. 
Since Aa p scales inversely with x to a steep power in the 
sub-Hill regime (sec [213), the contribution to the RMS 
speed is greatest from objects at small x (for reasonable 
variations of S TO with x). We take disk material to ex- 
tend to a minimum distance of |a; m inj = TvL-Rh (71 > 1) 
from the planet's orbit. Insertion of (|25() into (|39|l yields 



(«p,rnd) 



1/2 



• \ 1/2 

1 / £ m a£m \ R H 



(fi p At)V2 I Ml 



1 / X m a p m 



if 1< TZ < (vh/u) 1 / 2 - 

1/2 



evu, 



(fipAt)V2 ^ Ml 

if ft> (vn/u) 1 / 2 . 



(41) 

The RMS speed is maximized for 7t = 1. 

2.3.2. Super-Hill Velocities (u > Vu) 

Next we consider the noise generated by planetesimals 
with super- Hill random velocities (u > Vn). We refer to 
close encounters that change a p by the maximal amount, 
max|Aa p | ~ (m/M p )a p e (Equation |31|). as "maximal 
encounters." Maximal encounters, which occur at impact 
parameters b < GM p /u 2 , make an order-unity contribu- 
tion to the total super-Hill stochasticity. Non-maximal 
(more distant) encounters contribute to the total stochas- 
ticity through a Coulomb-like logarithm, as we show at 
the end of this sub-section. 

For a maximal encounter, a planetesimal must ap- 
proach within distance b < GM p /u 2 . Such encounters 
occur at a mean rate 



m V u J 



(42) 



where n ~ T, m il/(mu) is the number density of plan- 
etesimals, and we have assumed that planetesimal incli- 
nations and eccentricities are of the same order. Over 
a time interval At, the planet encounters on average 

NAt such planetesimals, each of which increases or de- 
creases dp by about max|Aa p |. Since planetesimals suf- 
fering maximal encounters have their orbits effectively 
randomized relative to each other, we may choose At 
to be any time interval longer than an encounter time 
At c ~ b/u < I /Sip (see related discussion in W2. 3.1() . 
Therefore the planet random walks with RMS velocity 
(averaged over time At) 



(^p,rnd) 



1/2 (ATAf) 1 / 2 max|Aa p 
At 



n p At 



1/2 



ia 2 m 



1/2 



M 2 



(43) 



grations (not shown). Since (d prnd 



and N cx b 2 , we have (d prnd 



What about the contribution from non-maximal en- 
counters? For a super-Hill encounter at impact parame- 
ter b, the specific impulse imparted to the planetesimal is 
~GM p /(bu). We suppose that |Aa p | is proportional to 
this specific impulse, so that |Aa p | oc 1/6. We have con- 
firmed this last proportionality by numerical orbit inte- 

i V2 cx (Ny/ 2 \Aa p \ 

1/2 

cx 6°, which implies 

that each octave in impact parameter contributes equally 
to the total stochasticity. In other words, our estimate 

for ( d 2 rnd ) in Equation l)43|l should be enhanced by 

a logarithmic factor of ln(6 max /6 m i n ), where 6 max and 
frmin ~ GM p /u 2 are maximum and minimum impact pa- 
rameters. We estimate 6 max ~ u/Q, the value for which 
the relative velocity of a super-Hill encounter is domi- 
nated by the planetesimal's random velocity rather than 
by the background shear. The logarithm is not large; for 
example, for e = 0.2, \n[(u/9)/(GM p /u 2 )] - 5. 

2.3.3. Summary 

We can neatly summarize Equations ((41(1 and 1(43(1 by 
defining the Hill eccentricity, 

e H = Rn/ap , (44) 

and parameterizing S m such that the disk contains mass 
M.AIp in planetesimals of mass m spread uniformly from 
a d /2 to 3a d /2: 



MMp 
27ra 2 



(45) 



where M is a dimensionless number of order unity. Then 

i /2 (Mm^ 1/2 



/•2 \ 1 /2 
\ a p,rnd/ 



CTZ- 4 (fipAt) -17 



V M p j a d 
if e < e K /II 2 



CK~ 2 (Q p At) 



-1/2 



[Mi 



1/2 



a d 



if en/71 2 < e < Tien; 



C (OpAt) 



1/2 



-i/2 ( Mm 

\ Mp J a d ( 
if e > 7?.eH 



(46) 

where we have introduced a constant coefficient C (the 
same for each case so that the function remains continu- 
ous across case boundaries). The coefficient C encapsu- 
lates all the factors of order unity that we have dropped 
in our derivations. By studying N-body simulation data 
pertaining to the case e > Tien as recorded in the lit- 
erature (Hahn & Malhotra 1999; Gomes et al. 2004), 
we estimate that C is possibly of the order of several. 
That C > 1 (but not 3> 1) is consonant with our hav- 
ing consistently underestimated the noise by neglecting 
(a) distant, non-maximal encounters in the super-Hill 



Stochastic Migration 



7 



regime ( H2. 3.2(1 . and (b) stochasticity introduced by ran- 
dom sampling of orbital elements of planctcsimals en- 
countering the planet f H2.3.1[l . We defer definitive cali- 
bration of C to future study, but retain the coefficient in 
our expressions below to assess the degree to which our 
quantitative estimates are uncertain. 

Planetesimals with e = TZe^i (|x m i n | = a p e) produce 
the maximum possible stochasticity: 

V2 „ 



:« rnd > 1/2 ~|(OpA*)- 1/2 ( 



Mm 



— e H «h, 



for e = 7£eH > en • 



(47) 



We will often adopt this case for illustration purposes 
below. 

/ v 1/2 

Since (a^ rnd \ cx (Mm) 1 ' , the stochasticity driven 

by planetesimals having a range of sizes is dominated by 
those objects (in some logarithmic size bin) having max- 
imal S m m. For common power-law size distributions, 
such objects will occupy the upper end of the distribu- 
tion. We will explicitly consider various possible size 
distributions in 



3. APPLICATION: MAXIMUM PLANETESIMAL SIZES 

Neptune is thought to have migrated outward by scat- 
tering planctcsimals during the late stages of planet for- 
mation (Fernandez & Ip 1984; Hahn & Malhotra 1999). 
As the planet migrated, it may have captured Kuiper belt 
objects into its exterior resonances (Malhotra 1995), giv- 
ing rise to the Resonant KBOs observed today (Chiang 
et al. 2003; Hahn & Malhotra 2005). If Neptune's mi- 
gration had been too stochastic, however, resonance cap- 
ture could not have occurred. A planetesimal disk having 
fixed surface mass density £ generates more stochasticity 
when composed of larger (fewer) planetesimals. There- 
fore, assuming that planetary migration and concomitant 
resonance capture correctly explain the origin of present- 
day Resonant KBOs, we can rule out size distributions 
that are too "top-heavy" during the era of migration. 

Stochasticity causes a planet to migrate both outward 
and inward. In H3.ll we provide background information 
regarding how resonance capture and retention depend 
on the sign of migration. In H3.2I we lay out general con- 
siderations for whether a stochastically migrating planet 
can retain particles in resonance. In H3.3I we derive and 
evaluate analytic, order-of-magnitude expressions for the 
maximum planetesimal size compatible with resonance 
retention, in the simple case when all planetesimals have 
the same size. In H3.4I we provide an analytic formula 
that details precisely how the resonance retention effi- 
ciency varies with average migration speed and planetes- 
imal size. These analytic results are quantitatively tested 
by numerical integrations in M3.5I Cases where planetes- 
imals exhibit a wide range of sizes are examined in H3.6I 

3.1. Migrating Outward and Inward 

As a planet migrates smoothly outward (away from the 
parent star) , it can capture planetesimals into its exterior 
mean-motion resonances. By contrast, a planet which 
migrates smoothly inward cannot capture planetesimals 
which are initially non-resonant into its exterior reso- 
nances (e.g., Peale 1986). But a planetesimal that starts 
in exterior resonance with an inwardly migrating planet 
can remain in resonance for a finite time. 




0.2 0.4 0.6 
Time (Myr) 



0.8 



Fig. 1. — Evolution of a planetesimal in external 3:2 resonance as 
the planet migrates smoothly inward. The planetesimal remains in 
resonance until its eccentricity reaches zero, at a semi-major axis 
of a = a min . 



Goldreich (1965) demonstrates how an outwardly mi- 
grating body adds angular momentum to a test particle 
in exterior resonance at just the right rate to keep the 
particle in resonance. Reversing the signs in his proof im- 
plies that an inwardly migrating body removes angular 
momentum from a particle in exterior resonance, pulling 
it inward while preserving the resonant lock. A planetes- 
imal's eccentricity decreases as it is pulled inward. The 
adiabatic invariant, 



TV = y/GM,a(py/l - e 2 - q) 



(48) 



which is preserved for migration timescales long com- 
pared to the synodic time (e.g., Murray-Clay & Chiang 
2005), implies that a planetesimal in p : q exterior reso- 
nance (p > q) cannot be pulled inward to a semi-major 
axis less than 



1 



GM* 



p-t 



«0 



o 



p-q 



the value for which e = 0. Here ao and eo are the ini- 
tial semi-major axis and eccentricity of the planetesimal, 
respectively. 

Thus an exterior particle follows an inwardly migrating 
planet in resonant loekstep until it cither reaches zero 
eccentricity (view Figure 4 of Peale 1986 or Figure 8.22 
of Murray & Dermott 1999 in reverse) or until it crosses 
the separatrix (view Figure 5 of Peale 1986 or Figure 8.23 
of Murray & Dermott 1999 in reverse), whichever comes 
first. We illustrate the former possibility in Figure ^and 
the latter possibility in Figure |3 using our own orbit 
integrations. The value of a m i n is annotated for reference. 
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Fig. 2. — Evolution of a planetesimal in external 3:2 resonance 
as the planet migrates smoothly inward. The planetesimal remains 
in resonance until it crosses the separatrix; by contrast to the evo- 
lution shown in Figure the particle's semi-major axis does not 
reach a mln and its eccentricity stays greater than zero. 



3.2. General Considerations for Resonance Retention 

The random component of the planet's migration is a 
form of Brownian motion. The planet encounters a large 
number of small planetesimals, each of which causes the 
planet's semi- major axis a p to randomly step a small 
distance. Each random change in a p produces a corre- 
sponding change in the semi-major axis of exact reso- 
nance, 7 but no corresponding change in the actual semi- 
major axis of a resonant particle. The semi-major axis 
of the particle does not respond because a p changes ran- 
domly every encounter time At e , which is much shorter 
than the resonant libration period. Only changes in a p 
that are coherent over timescales longer than the libra- 
tion period produce an adiabatic response in the parti- 
cle's semi-major axis. We have verified these assertions 
by numerical orbit integrations (not shown). 

The no-response condition implies that as a p random 
walks, the difference between the semi-major axes of ex- 
act resonance and of a resonant particle random walks 
correspondingly. In other words, a resonant particle's 
libration amplitude random walks. The sign and mag- 
nitude of each step in the libration amplitude's random 
walk depend on the phase of libration when the step is 
taken. Since at any given time an ensemble of resonant 
particles are distributed over the full range of phases, a 
single random-walk history for the planet generates an 
ensemble of different random-walk histories for the par- 
ticles. 

When the libration amplitude of a resonant particle 

7 A particle in exact resonance has zero libration amplitude, by 
definition. 



random walks past its maximum allowed value, the par- 
ticle escapes resonance. The maximum libration ampli- 
tude (full width) as measured in semi-major axis is 

& P ,iib = 2C libap ( ^I2l) 1/2 , (50) 



where e rGS is the eccentricity of the resonant object (not 
to be confused with the planetesimals generating the bulk 
of the noise), Cub ~ 4^// 3 i/3 is a constant (see Murray 
and Dermott 1999 for /31), and we have restricted con- 
sideration to first-order (p — q = 1) resonances. For the 
3:2 exterior resonance with Neptune, Cub ~ 3.64. Note 
that in contrast to the usual definition of maximum libra- 
tion width, (5a p jib refers not to the particle's semi-major 
axis, but rather to the planet's. The meaning of <5a p jib is 
as follows. Take a particle in exact resonance. By defini- 
tion, such a particle has zero libration amplitude. Then 
the planet's semi-major axis can change instantaneously 
by at most 5a Pl nh/2 and the particle will still remain in 
resonance (but with finite libration amplitude). 

Equation (|50|l derives from the pendulum model of res- 
onance, which is known to be inaccurate at large e res for 
some resonances. Malhotra (1996) finds numerically that 
for e rcs = 0.1-0.4, ^a p j;b for the 3:2 resonance is insen- 
sitive to e rcs , whereas the pendulum model predicts that 
^Opjib doubles over this range. We nevertheless employ 
Equation l|50|) to estimate the maximum libration width, 
since it is simple, analytic, and introduces errors less than 
of order unity in our numerical evaluations below. The 
qualitative physics described in this paper does not de- 
pend on the accuracy to which we estimate <5a Pj iib- 

Consider a planet which migrates outward on aver- 
age. When the random component of the planet's mi- 
gration is added to the average component, a planet 
can migrate either outward or inward at any moment. 
Call iSmd = f <jp,rnd dt the running sum of the random 
changes in a p . The probability Pkeep that a given par- 
ticle is retained in resonance over some duration of mi- 
gration equals the probability that ISVndl remains less 
than the maximum libration half- width <5a p jib/2 during 
that time. A particle that escapes resonance by being 
dropped behind the resonance (tSVnd — d^^pjib/z) is, 
practically speaking, permanently lost. The planet can- 
not recapture the particle by smoothly migrating inward 
(see Jl . The random component of the planet's mi- 
gration can cause the planetesimal to be recaptured, but 
a recaptured particle lies on a trajectory near the separa- 
trix and quickly re-escapes in practice. Once the average 
(outward) component of the planet's migration carries 
the resonance well past the particle, the particle cannot 
be recaptured even if SVnd random walks back to zero; in 
other words, the particle has been permanently left be- 
hind. A particle that escapes by being dropped in front 
of the resonance (Smd = — &Jp,iit>/2) is also lost more of- 
ten than not. Such a particle can be recaptured when the 
planet resumes migrating outward. Nevertheless, upon 
its recapture onto a trajectory near the separatrix, the 
particle can librate back to smaller semi-major axes and 
be expelled behind the resonance permanently. 

3.3. Order-of-Magnitude Planetesimal Sizes 

Armed with the considerations of ^ 13.21 we are now 
ready to derive analytic, order-of-magnitude expressions 
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for the maximum planetesimal sizes compatible with res- 
onance retention, for the simple case when the disk is 
composed of objects of a single size. The assumption of 
a single size is relaxed in H3.6I 

Say the planet takes time T to migrate at speed a Pj avg 
from its initial to its final semi-major axis. Over this 
time, a p random walks an expected distance of 

,1/2, 



' (<2p,rnd) T 



_ 4 / Mm 



C7Z 



Mm 



1/2 



— e H u H 



T 



1/2 



if e < e H /7Z 2 



1/2 



■ev H 



1/2 



T 

a d ' (51) 

if en/11 2 <e< fte H ; 



(Mi 



1/2 



vn 

ad e 



T 



1/2 



if e > TZen, 



l p,rnd 



1/2 



where we have set At = T in evaluating 

If °a p ,T < Sa Pj nb/2, the planet can keep a large fraction 
of planetesimals in resonance. That is, most particles 
are retained in resonance when the disk mass comprises 
planetesimals of mass 

\ 2 i 

\ -L ^res * 

ilpT e H 



K s C 2 h 
C 2 M 



if e < e H /7Z 2 ; 



m < m crit ~ < 



^C 2 h 
C 2 M 



r 2 
°iib 



s eH 



rtpT 



if en/71 2 <e< Tien; 



C 2 M 



«d 



I e e 2 



flpT 



if e > IZen. 

(52) 

Equation 1(521) can be equivalently interpreted as an up- 
per limit on T for planetesimals of given mass to. For a 
fixed degree of noise, resonant objects are more difficult 
to retain if the average migration is slow. 

We evaluate (|52f) to estimate the maximum planetes- 
imal radius, s — (3m/47rp) 1 / 3 , compatible with reso- 
nant capture of KBOs by Neptune. For an internal 
density p = 2g/cm 3 , M p = M N = 17M e , e H = 0.03, 
dp = a d = 26.6 AU, e ros = 0.25, M = 2 (so that 
S m = 0.2 gem -2 ), and T — 3 x 10 7 yr, resonant cap- 
ture and retention require 

700ft 8 / 3 C- 2 / 3 , ife<e H /ft 2 ; 

70e- 2 / 3 ^ 4 / 3 C~ 2 / 3 , if e H /K 2 < e < Ke H ; 

7000e 2 / 3 C- 2 / 3 , ife>"Ke H . 

(53) 

For example, if e = 0.1 and 71=1, then line 3 of (|53|l ob- 
tains and s cr ;t = 1500C -2 / 3 km. Maximum stochasticity 



km 



< 



^crit 

km 



results when TZ = 1 and e < en (see also Equation |47j): 
either of lines 1 or 2 then yield s cr ;t ~ 700C _2 ' 3 km. 
These size estimates decrease by about 20% when cor- 
rected to reflect the fact that the width of the 3:2 res- 
onance is somewhat smaller than the pendulum model 
implies (see the discussion following Equation 150(1. 

3.4. Analytical Formula for the Retention Fraction 

As defined in i !3.2l Pk eep is the resonance retention frac- 
tion, or the probability that a typical resonant particle 
is retained in resonance over some duration of migration. 
We calculate Pkeep by modelling the random component 
of the planet's migration as a diffusive continuum pro- 
cess. In the limit that the planet encounters a large 
number N 3> 1 of planetesimals, the Poisson distribu- 
tion (Equation |37| 1 is well-approximated by a Gaussian 

distribution with mean N and variance N. Thus, over 
■ -l 

a time interval At 3> N , the random displacement of 
the planet, ASVnd = Aa p (N — TV), has the probability 
density distribution 

f(AS md ,At) = -=L=exp(-(A5 rnd ) 2 /(2 J DAt)) , 
V zirDAt 

(54) 

where D = (Aa p ) 2 N is the diffusion coefficient and we 
recall that Aa p is the change in a p due to an encounter 
with a single planetesimal. The evolution of a P) rnd with t 
is continuous and the distribution / is independent over 
any two non-overlapping intervals At (the random walk 
has no memory). In other words, AS^d evolves as a 
Wiener process, or equivalently according to the rules of 
Brownian motion (e.g., Grimmett & Stirzaker 2001a). 
From Equation 1(54(1 . it follows that over time T, the 
probability that | ASVnd] does not exceed <5a Pj nb/2 equals 



Pkeep — 



n=l 



4 , /nn 
— sin 3 — 
nir V 2 



\ p -\ n T 



(55) 



where A„ = (nn) 2 D / (2Sa 2 lih ) (see Appendix IaI for a 
derivation) . 

Suppose migration occurs in a disk of planetesimals 
having a single size s. Figure [3] displays Pkeep as a func- 
tion of s and of exponential migration timescale r defined 
according to 



5 (*) 



..i)e 



-t/r 



(56) 



a P ,f — (a>p,{ 

where a Pji and a p [ are the planet's initial and final aver- 
age semi- major axes, respectively. In Equation 1)55(1 ■ we 
take T = 2.6r, and evaluate remaining quantities for the 
case of maximum stochasticity: e < en and TZ = 1. Then 

Aa p = Aopjci (Equation j25]) and TV = 2Y, rn flR 2 1 /m 
(Equation |36( , with a factor of 2 inserted to account for 
disk material both inside and outside the planet's orbit). 



As in $121 we take p = 2g/cm 3 , M p 



17 Mr, 



and M 



en = v. 06, a p = a d = Zb.b AU, Kn = ena p , 
To evaluate <5a Pj i;b, we take e rcs = 0.25 for a particle in 
3:2 resonance. Figure |3 describes how for a given size 
s, the retention fraction decreases with increasing t; the 
longer the duration of migration, the more chance a par- 
ticle has of being jostled out of resonance. For r = 10 
Myr, planetesimals must have sizes s < 500 km for the 
retention fraction to remain greater than 1/2. These re- 
sults confirm and refine our order-of-magnitude estimates 
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made in H3.3I Similar results were obtained for the 2:1 
resonance. 

The continuum limit is valid as long as the expectation 
value of the time required for a resonant particle to es- 

— _1 2 

cape, (£ csca pc) ~ N [((5a P! iib/2)/Aa p ] , greatly exceeds 

the time for the planet to encounter one planetesimal, 
■ -l 

N . This criterion is satisfied for the full range of pa- 
rameters adopted in Figure |3| 

3.5. Numerical Results for the Retention Fraction 

To explore how a stochastically migrating planet cap- 
tures and retains test particles into its exterior reso- 
nances, and to test the analytic considerations of S H3.2I - 
13.41 we perform a series of numerical integrations. We 
focus as before on the 3:2 (Plutino) resonance with Nep- 
tune. 

Following Murray-Clay & Chiang (2005, hereafter 
MC05), we employ a series expansion for the time- 
dependent Hamiltonian, 

1/2 

(2r 



H = - 



(GM g 
2(3r 
GM n 



-AO 2 

Hh - 



GMr. 



he 2 + /3iecos< 



■AO 



(57) 



a p (t) 

where a — a p /a w 0.76, the /Vs are given in Murray & 
Dermott (1999), and Af (Equation |48|) is a constant of 
the motion determined by initial conditions. The reso- 
nance angle, 

4> = 3A rcs — 2A P — n7 rcs , (58) 

is defined by the mean longitude A rcs and longitude of 
periastron vj Tes of the resonant particle, and the mean 
longitude A p of the planet. The resonance angle librates 
about 7r for particles in resonance. The momentum con- 
jugate to <j> is r. We integrate the equations of motion, 

■ dH ■ dH 

+ = W> T = ~W m 

using the Bulirsch-Stoer algorithm (Press et al. 1992) for 
fixed a and fa's. 

The Hamiltonian in Equation l|57|l faithfully repro- 
duces the main features of the resonance potential; see 
Beauge (1994, his Figures 12a and 12c) for a direct com- 
parison between such a truncated Hamiltonian and the 
exact Hamiltonian, averaged over the synodic period (see 
also that paper and Murray-Clay & Chiang 2005 for a 
discussion of the pitfalls of keeping one too many a term 
in the expansion). Of course, even the exact Hamilto- 
nian, because it is time-averaged and neglects chaotic 
zones, is inaccurate with regards to details such as the 
libration width, but these inaccuracies are slight; see the 
discussions following Equations H50(l and 1)53(1 . 

To compute a p (i), we specify separately the average 
and random components of the migration velocity, a P ,avg 
and a P! rnd- For d p avK , we adopt the prescription (equiv- 
alent to Equation 



1 



-^avg 



(a 



-t/r 



(60) 



where a Pi j and a Pi f are the planet's initial and final av- 
erage semi-major axes, respectively, and r is a time con- 
stant. To compute d Pirnc i, we divide the integration into 



time intervals of length l/fl p . The only requirement for 
the time interval is that it be less than the libration pe- 
riod tub ~ 400/f2 p (see il3.2p . Over each interval, we 
randomly generate 

a p , rnd = n p Aa p {N n - M^) . (61) 

We focus on the case of maximum stochasticity, so that 

Aa p = Aa P: nc i (Equation [21]) and N — 2£ m i? 2 i f2 p /m 
(Equation |36j with 1Z = 1 and an extra factor of 2 in- 
serted to account for disk material on both sides of the 
planet's orbit). We assume that the entirety of the disk 
mass is in planetesimals of a single mass to. Each Nq 
is a random deviate drawn from a Poisson distribution 
having mean Nil^ 1 . 

Figure^displays the sample evolution of a test particle 
driven into 3:2 resonance by a stochastically migrating 
planet. For this integration, M p — Mn, a p .i = 23.1 AU, 
ap, f = 30.1 AU, t = 10 7 yr, K = 1, M = 2 (so that 
E m = 0.2gcm~ 2 ), and s = 150km (to = 3 x 10 22 g). 
This choice for s is sufficiently small that the particle is 
successfully captured and retained by the planet. 

Contrast Figure 01 with Figure in which all model 
parameters are the same except for a larger s = 700 km. 
In this case the planet eventually loses the test particle 
because the migration is too noisy. 

Figure El displays the fraction of particles caught and 
kept in resonance as a function of r and s. For each data 
point in Figure we follow the evolution of 200 parti- 
cles initialized with eccentricities of approximately 0.01 
and semi-major axes that lie outside the initial position 
of resonance by about 1 AU. 8 Figure is the numerical 
counterpart of Figure the agreement between the two 
is excellent and validates our analytic considerations. If 
r = 10 7 yr (consistent with findings by MC05), then the 
capture fraction rises above 0.5 for s < 500 km. Since 
these results pertain to the case {1Z = 1, e < en} which 
yields the largest amount of noise for given S m and s, 
we conclude that s ~ 500 km is the lowest, and thus the 
most conservative, estimate we can make for the maxi- 
mum planetesimal size compatible with resonant capture 
of KBOs by a migrating Neptune, assuming that the en- 
tire disk is composed of planetesimals of a single size 
(this assumption is relaxed in the next section). In other 
words, if Neptune's migration were driven by planetesi- 
mals all having s <ti 500 km, stochasticity would not have 
impeded the trapping of Resonant KBOs. Of course, 
our numerical estimate of 500 km is uncertain insofar as 
we have not kept track of order-unity constants in our 
derivations. We suspect a more careful analysis will re- 
vise our size estimate downwards by a factor of a few (see 
the discussion of C in M2.3.3JI . 

3.6. Planetesimal Size Distributions 

Actual disks comprise planetesimals with a range of 
sizes. From Equation H46|l. the stochasticity in the 
planet's migration is dominated by those planetesimals 
having maximal H m m. What was the distribution of 
sizes during the era of Neptune's migration? A possi- 
ble answer is provided by the coagulation simulations 

8 The particles do not all have the same initial eccentricities and 
semi-major axes. This is because they occupy the same Hamilto- 
nian level curve; see section 3.5 of MC05. 
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Fig. 3. — Fraction of particles re tain ed in external 3:2 resonance by a stochastically migrating planet as a function of migration timescale, 
calculated according to Equation 1551 , The entire disk mass is assumed to be in planetesimals of a single size s, and a range of choices for 
s are shown. The diffusivity D is evaluated at its maximum value, appropriate for the case e < ejj and 1Z = 1. We set Aa p = Aa P;nc i 

(Equation |2!>]), N = 2£ m QR 2 J /m, (Equation p = 2g/cm 3 , M p = M N = 17M e , e H = 0.03, a p = a d = 26.6 AU, R H = e H a p , 

M = 2, e res = 0.25, and T = 2.6t. Planetesimals having sizes smaller than ~200 km produce so little noise in the planet's migration that 
no object is lost from the 3:2 resonance. Compare this Figure with its numerical counterpart, Figure [5] In calculating Pkeepi we assume 

C = 1; probab ly C is of order several, in which case the sizes indicated in the Figure should be revised downward by a factor of a few (C 2 / 3 ; 
see Equation 1531 1. 



of Kenyon & Luu (1999, hereafter KL99). The left- 
hand panel of their Figure 8 portrays the evolution of 
the size distribution, starting with a disk of seed bodies 
having sizes up to 100 m and a total surface density of 
£ = 0.2gcm~ 2 . After t = 11 Myr, the size bin for which 
S m m is maximal is centered at s ~ 4 km; for this bin at 
that time, E m = 10~ 3 gcm -2 (evaluated within a loga- 
rithmic size interval 0.3 dex wide). After t = 37 Myr, 
the planetesimals generating the most stochasticity have 
s ~ 750 km and S m = 2 x 10~ 3 gcm~ 2 . Note that at 
t = 37 Myr, the stochasticity is dominated by the largest 
planetesimals formed, but they do not contain the bulk 
of the total disk mass; the lion's share of the mass is 
instead sequestered into km-sized objects. 

In Figure [7] we plot the resonance retention fraction 
Pkccp (Equation [55]) for the KL99 size distribution at 
t = 11 and 37 Myr, using the values of s(m) and E m 
cited above. The remaining parameters that enter into 
Pkeep are chosen to be the same as those employed for 
Figure|3| i.e., we adopt the case of maximum stochastic- 
ity. Evidently, Pkccp = 1 for the KL99 size distributions; 
stochasticity is negligible. 

For comparison, we also plot in Figure the retention 
fraction for pure power-law size distributions: drj/ds oc 
s~ q , where drj is the differential number of planetesimals 
having sizes between s and s + ds. Since S m m oc s 7 ~ q , 
stochasticity is dominated by the upper end of the size 



distribution for q < 7. We fix the maximal radius to be 
that of Pluto (s U ppcr = 1200km), set the total surface 
density £ = 0.2 gem -2 , and calculate Pkccp for three 
choices of q = 3.5, 4, and 4.5. For q > 4, the lower 
limit of the size distribution significantly influences the 
normalization of drj/ds] for q = 4 and 4.5, we exper- 
iment with two choices for the minimum planetesimal 
radius, si owor = 1km and lm. We equate E m with the 
integrated surface density between s upP er/2 and s upP er- 
According to Figure steep size distributions q > 4 
are characterized by order-unity retention efficiencies. In 
contrast, shallow size distributions q < 4 for which the 
bulk of the mass is concentrated towards s uppcr can in- 
troduce significant stochasticity. 

4. CONCLUDING REMARKS 

We summarize our findings in H4.1l and discuss quanti- 
tatively some remaining issues in H4.2I 

4.1. Summary 

Newly formed planets likely occupy remnant planetes- 
imal disks. Planets migrate as they exchange energy and 
angular momentum with planetesimals. Driven by dis- 
crete scattering events, migration is stochastic. 

In our solar system, Neptune may have migrated out- 
ward by several AU and thereby captured the many 
Kuiper belt objects (KBOs) found today in mean-motion 
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Time (Myr) 

Fig. 4. — Evolution of a particle caught into 3:2 resonance with a 
stochastically migrating planet. Stochasticity is driven by a disk of 
surface density E m = 0.2 gem -2 , all in planetesimals having sizes 
s = 150 km and sub-Hill random velocities. The random walk in 
the planet's semi-major axis causes the libration amplitude of the 
resonant particle to undergo a corresponding random walk. The 
noise in this example is too mild to prevent the planet from both 
capturing and retaining the particle in resonance. 
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Fig. 5. — Evolution of a particle caught into, but even- 
tually lost from, 3:2 resonance with a stochastically migrating 
planet. Stochasticity is driven by a disk of surface density S m = 
0.2 gem -2 , all in planetesimals having sizes s = 700km and sub- 
Hill random velocities. The particle is expelled from resonance 
having had its eccentricity raised to 0.2 during its time in resonant 
lock. 



resonance with the planet. While resonance capture is 
efficient when migration is smooth, a longstanding is- 
sue has been whether Neptune's actual migration was 
too noisy to permit capture. Our work addresses — and 
dispels — this concern by supplying a first-principles the- 
ory for how a planet's semi-major axis fluctuates in re- 
sponse to intrinsic granularity in the gravitational poten- 
tial. We apply our theory to identify the environmental 
conditions under which resonance capture remains vi- 
able. 

Stochasticity results from random variations in the 
numbers and orbital properties of planetesimals encoun- 
tering the planet. The degree of stochasticity (as mea- 
sured, say, by c 0pi T, the typical distance that the planet's 
semi-major axis random walks away from its average 
value) depends on how planctesimal semi-major axes a 
and random velocities u are distributed. We have param- 
eterized a by its difference from the planet's semi-major 
axis: x = a — a p = TZRu, where i?H is the Hill sphere 
radius and 1Z > 1. In the case of high dispersion when 



u > IZvft (where vjj = ^ p ^?h is the Hill velocity and 
Sip is the planet's orbital angular velocity), planetesimal 
orbits cross that of the planet. Stochasticity increases 
with decreasing u in the high-dispersion case because 
the cross-section for strong scatterings increases steeply 
with decreasing velocity dispersion (as 1/w 4 )- In the 
intermediate-dispersion case when vh/1Z 2 < u < TZva, 
planetesimal and planet orbits do not cross, and stochas- 
ticity decreases with decreasing u. In the low-dispersion 
case when u < vn/TZ 2 , the amount of stochasticity is 
insensitive to u. 

The values of u and 1Z which actually characterize disks 
are unknown. The random velocity u, for example, is 
expected to be set by a balance between excitation by 
gravitational scatterings and damping by inelastic colli- 
sions between planetesimals and/or gas drag. Damping 
depends, in turn, on the size distribution of planetesi- 
mals. These considerations are often absent from current 
N-body simulations of planetary migration in planetesi- 
mal disks. Despite such uncertainty, we can still identify 
the circumstances under which stochasticity is maximal. 
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Fig. 6. — Fraction of particles caught into, and retained within, external 3:2 resonance by a stochastically migrating planet. For every t 
and s, we numerically integrate the trajectories of 200 test particles with initial eccentricities of ~0.01 and semi-major axes that lie 1 AU 
outside of nominal resonance. These particles respond to the time-averaged potential of a Neptune-mass planet which migrates outward 
from 23.1 AU to 30.1 AU within a disk of fixed surface density S m = 0.2 gem -2 in planetesimals of a single size s. The planetesimals 
have sub-Hill random velocities and semi-major axes that lie within TZ. = 1 Hill radius of the planet's; these choices maximize the amount 
of stochasticity in the planet's migration. Compare this Figure with its analytic counterpart, Figure [3] the agreement is excellent. The 
solid curve labelled "Smooth" corresponds to the case when all noise is eliminated from the planet's migration. Planetesimals having sizes 
smaller than ~200 km yield an essentially smooth migration. For r < 10 5 yr, capture is not possible even if migration were smooth, since 
the migration is too fast to be adiabatic. These results are calculated for C = 1; probably C is of order several and so the sizes indicated in 
the Figure should be revised downward by a factor of a few (C 2 / 3 ; see Equation I53I - ). 



Maximum stochasticity obtains when 1Z ~ 1 and u < «h, 
that is, when planetesimals have semi-major axes within 
a Hill radius of the planet's and when their velocity dis- 
persion is no greater than the Hill velocity. 

A stochastically migrating planet cannot retain objects 
in a given resonance if the planet's semi-major axis ran- 
dom walks away from its average value by a distance 
greater than the maximum libration width of the reso- 
nance. This simple criterion is validated by numerical 
experiments and enables analytic calculation of the res- 
onance retention efficiency as a function of disk param- 
eters. A disk of given surface density generates more 
noise when composed of fewer, larger planetesimals. In 
the context of Neptune's migration, we estimate that if 
the bulk of the minimum-mass disk resided in bodies 
having sizes smaller than 0(100) km and if the fraction 
of the disk mass in larger bodies was not too large (< 
a few percent for planetesimals having sizes of 1000 km, 
for example), then the retention efficiency of Neptune's 
first-order resonances would have been of order unity 
(> 0.1). Such order-unity efficiencies seem required by 
observations, which prima facie place 122/474 ss 26% of 
well-observed KBOs (excluding Centaurs) inside mean- 
motion resonances (Chiang et al. 2006). Drawing conclu- 
sions based on a comparison between this observed per- 
centage and our theoretical retention percentage Pkeep 



is a task fraught with caveats — a more fair comparison 
would require, e.g., disentangling the observational bias 
against discovering Resonant vs. non-Resonant objects; 
account of the attrition of the Resonant population due 
to weak chaos over the four-billion-year age of the so- 
lar system; and knowledge of the initial eccentricity and 
semi-major axis distributions of objects prior to reso- 
nance sweeping, as these distributions impact capture 
probabilities in different ways for different resonances 
(Chiang et al. 2003; Hahn & Malhotra 2005; Chiang et 
al. 2006). But each of these caveats alters the relevant 
percentages only by factors of a few, and when combined, 
their effects tend to cancel. Therefore we feel comfort- 
able in our assessment that -Pkeep must have been of or- 
der unity to explain the current Resonant population. In 
that case, O(100 km) is a conservative estimate for the 
maximum allowed size of planetesimals comprising the 
bulk of the disk mass, derived for the case of maximum 
stochasticity. 

How does an upper limit of 0(100) km compare with 
the actual size distribution of the planetesimal disk? 
While today's Kuiper belt places most of its mass in 
objects having sizes of ~100 km, this total mass is 
tiny — only ^{).\M§ (Bernstein et al. 2004; see Chiang et 
al. 2006 for a synopsis). The current belt is therefore 2-3 
orders of magnitude too low in mass to have driven Nep- 
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Fig. 7. — Fraction of particles retained in external 3:2 resonance by a stocha stic ally migrating planet for various planetesimal size 
distributions. The retention efficiency is calculated analytically using Equation 1551 . with parameters the same as those for Figure 151 
except for S m X m; that parameter is evaluated at its maximum value within a logarithmic size bin spanning a factor of 2 for a given size 
distribution. The size distributions considered include two from Kenyon & Luu (1999; their Figure 8), evaluated at times t = 11 Myr and 
37 Myr; and five different power-law distributions, each characterized by a total integrated surface density S = 0.2 gem -2 , an upper size 
limit Supper = 1200 km, a differential size index q (such that drj/ds oc s~ 9 ), and a lower size limit Slower as indicated (the curve for q = 3.5 
is insensitive to si OWO r since the bulk of the mass is concentrated towards s upP er). The three curves for the size distributions of KL99 and 



for {q = 4.5, s\ 



1 m} overlap at P\ 



keep 



tune's migration. The current size distribution is such 
that bodies having radii > 40 km are collisionless over the 
age of the solar system and might therefore represent a 
direct remnant, unadulterated by erosive collisions, of the 
planetesimal disk during the era of migration (Pan & Sari 
2005). If so, the bulk of the primordial disk mass must 
have resided in bodies having sizes < 40 km. Theoreti- 
cal calculations of the coagulation history of the Kuiper 
belt are so far consistent with this expectation. Kenyon 
& Luu (1999) find, for their primordial trans-Neptunian 
disk of 10M®, that 99% of the mass failed to coagu- 
late into bodies larger than 0(1) km, because the forma- 
tion of several Pluto-sized objects (comprising ^0.1% of 
the total mass) excited velocity dispersions so much that 
planetesimal collisions became destructive rather than 
agglomerative. The average-mass planetesimals in their 
simulation have sizes 0(1) km, much smaller than even 
our most conservative estimate of the maximum allowed 
size of 0(100) km. 

For a given size distribution of planetesimals, most 
stochasticity is produced by the size bin having maxi- 
mal rj in , which need not be the size bin containing the 
majority of the mass. Here, rj and m are the number 
of planetesimals and the mass of an individual planetes- 
imal in a logarithmic size bin. For power-law size dis- 
tributions drj/ds oc s~ q such that q < 7, stochasticity is 
dominated by the largest planetesimals. For disks having 



as much mass as the minimum-mass disk of solids and 
whose largest members are Pluto-sized, size distributions 
with q > 4 enjoy order-unity efficiencies for resonance re- 
tention. The size distributions of Kenyon & Luu (1999) 
resemble q = 4 power laws, but with a large overabun- 
dance of planetesimals having sizes of 0(1) km. This 
sequestration of mass dramatically reduces the stochas- 
ticity generated by the largest bodies, which have sizes 
of 0(1000) km. 

We conclude that Neptune's Brownian motion did not 
impede in any substantive way the planet's capture and 
retention of Resonant KBOs. 

4.2. Extensions 

4.2.1. Single Kick to Planet 

Our focus thus far has been on the regime in which 
many stochastic kicks to the planet are required for res- 
onant particles to escape. Of course, a single kick from a 
planetesimal having sufficiently large mass mi could flush 
particles from resonance. To estimate mi , we equate the 
change in the planet's semi-major axis from a single en- 
counter, Acip, to the maximum half- width of the reso- 
nance, 5a p ,iib/2 (see Equation [SUj and related discus- 
sion). In the likely event that the perturber's eccentric- 
ity e is of order unity, then max(Aa p ) ~ (rai/M p )o p e 
(Equation j3J) and therefore m\ > 0.6(0.5/e)M® for 
Plutinos to escape resonance. Our estimate for m\ agrees 
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with that of Malhotra (1993). Why such enormous per- 
turbcrs have not been observed today is unclear and 
casts doubt on their existence (Morbidelli, Jacob, & Petit 
2002). If such an Earth-mass planetesimal were present 
over the duration T of Neptune's migration, then the 
likelihood of a resonance-destabilizing encounter would 

be Pi ~ NT ~ 10 _2 (0.5/e) 4 , where the encounter rate 
N is given by Equation l|42|l with S m ~ mi/(27TOp), and 
we have set T - 2.7 x 10 7 yr. 

4.2.2. Kicks to Resonant Planetesimals 

Finally, we have ignored in this work how disk planetes- 
imals directly perturb the semi-major axis of a resonant 
particle. This neglect does not significantly alter our con- 
clusions. Take the resonant planetesimal to resemble a 
typical Resonant KBO observed today, having size s rcs ~ 
100 km. Then its Hill velocity is en.res^a ~ 10 -4 fla. The 
relative velocity between the resonant planetesimal and 
an ambient, perturbing planetesimal greatly exceeds this 
Hill velocity, if only because migration in resonant lock 
quickly raises the eccentricity of the resonant planetes- 

9 The probability that a planetesimal will be ejected from reso- 
nance by a planetesimal of comparable mass in a single encounter 



imal above en. res- Equation i|43[l . appropriate for the 
super-Hill regime, implies that (<ip rn d) oc M — the 
RMS random velocity does not depend on the mass of the 
object being perturbed! Therefore when both the reso- 
nant planetesimal and the planet are scattering planetes- 
imals in the super-Hill regime, their random walks are 
comparable in vigor. The conservative limit of 0(100) 
km on the planetesimal size which allows resonance re- 
tention is derived, by contrast, for the sub-Hill, maxi- 
mum stochasticity regime, and is therefore little affected 
by these considerations. 9 
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APPENDIX 

ABSORPTION PROBABILITY FOR BROWNIAN MOTION WITH A DOUBLE BOUNDARY 

Here we derive Equation l|55|) , the probability that a particle experiencing Brownian motion between two absorbing 
boundaries has not been absorbed by time t (e.g., Grimmett and Stirzaker 2001b). Consider Brownian motion along 
a path x(t) with x(0) = and absorbing boundaries at x = ±6, b > 0. The probability density distribution f(x,t) 
satisfies the diffusion equation 



9f = l n ^f 

dt 2 dx 2 ' 

where D is the diffusion coefficient. The absorbing boundaries generate the boundary conditions 

/(±M) = o 

for all time, 10 and the initial condition is 

f(x,0) = 5(x) . 

To solve for f(x,t), we expand / in a Fourier series, keeping only terms that satisfy (|A2|) : 

f(x,t)=J2 k n(t) sinf 
n=l 

Plugging (fA"4)) into l[Alj) . we find 

k n (t) 

where the c n 's are constants and 

nV 

The c ra 's must satisfy <|A3ll . From Fourier analysis at time t — 0, we find 



2b 



-A„t 



1 

2b Jo 



5(y-b)-6(y-3b)] 



1 /mr\ 

r in h-J 



(Al) 

(A2) 
(A3) 

(A4) 

(A5) 
(A6) 

(A7) 
(A8) 



10 Equation IA2I holds as long as D is non-zero. Particles near the boundary are carried across by fluctuations too quickly to maintain 
a non-zero density / at x = ±fe (see Grimmett and Stirzaker 2001a for a proof). 
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The probability that the walker has not yet crossed either of the absorbing boundaries at time t is 

rb 



(*)=/ f(x,t)dx (A9) 
J-b 



b 

b oo 



E^-^sin(^±^)dx (A10) 



'n=l 



^± sin 3(!^) e -(-)^/(8 & 2 ). (AU) 
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estimated to be of order several 
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